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The mean-field dynamics of a Bose-Einstein condensate is studied in presence of a microscopic 
trapping potential from which the condensate can escape via tunneling through finite barriers. We 
show that the method of complex scaling can be used to obtain a quantitative description of this 
decay process. A real-time propagation approach that is applied to the complex-scaled Gross- 
Pitaevskii equation allows us to calculate the chemical potentials and lifetimes of the metastably 
trapped Bose-Einstein condensate. The method is applied to a one-dimensional harmonic confine- 
ment potential combined with a Gaussian envelope, for which we compute the lowest symmetric 
and antisymmetric quasibound states of the condensate. A comparison with alternative approaches 
using absorbing boundary conditions as well as complex absorbing potentials shows good agreement. 
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I. INTRODUCTION 

With the advent of optical lattices and "atom chips" 
PJ, it became possible to probe the transport properties 
of a Bose-Einstein condensate in the mesoscopic regime. 
Experiments in this context include the observation of 
Bloch oscillations 0, El- the guided and free propaga- 
tion of condensates through waveguide structures 0, || , 
the transport of condensates with "optical tweezers" [6j , 
as well as the realization of Josephson junctions Q and 
matter-wave interferometry [8j, to mention just a few ex- 
amples. Those experiments typically involve rather small 
trapping potentials, with length scales that can be of the 
order of a few microns. In such microscopic geometries, 
decay mechanisms of the condensate become a relevant 
issue. On the one hand, the condensed state is, at finite 
atom densities, subject to depletion, which is caused by 
the interaction with the thermal cloud and by three-body 
collisions. On the other hand, the condensate can escape 
from the trapping potential by tunneling though its bar- 
riers, if the chemical potential of the condensed atoms 
exceeds the background potential in the free space out- 
side the trap. In that case, the self-consistent mean-field 
state of the condensate is no longer bound, but rather 
corresponds to a metastable "resonance" state, in a sim- 
ilar way as, e.g., doubly excited electronic states in the 
helium atom. 

From the theoretical point of view, the problem of 
metastable states of Bose-Einstein condensates in such 
"open" trapping potentials was first approached by Moi- 
seyev et al. in Ref. 9]. In this work, the mean-field dy- 
namics of a condensate was studied in presence of an 
isotropic harmonic confinement that is combined with 
a Gaussian envelope. Decaying states of the conden- 
sates were identified with self-consistent solutions of the 
stationary Gross-Pitaevskii equation at complex-valued 
chemical potentials, where complex absorbing potentials 
at the boundaries of the numerical grid were used to ac- 
count for the escape of the condensate. As a result, a 
cross-over from a decaying (resonance) state to a bound 
state of the condensate was found for a finite attractive 



interaction between the atoms. 

The calculations of Ref. [i| were satisfactory from the 
quantitative point of view, but raised an important open 
question: How can "resonances", i.e. stationary states 
that describe the escape of population from an open con- 
finement potential, be formally introduced in the frame- 
work of the nonlinear Gross-Pitaevskii equation? For lin- 
ear systems, it is well known that this task is most con- 
veniently accomplished by applying the method of "com- 
plex scaling" (or "complex rotation") [Tol fill [T^. 
This technique essentially amounts to the complex di- 
lations r i — ► re 10 and — iV i— > — iVe~ l8 of the position 
and momentum operators in the Hamiltonian that de- 
scribes the quantum system under study. This transfor- 
mation leads to a nonhcrmitean Hamiltonian with a com- 
plex eigenvalue spectrum the continuous part of which is 
rotated to the lower half of the complex energy plane. 
Resonances, i.e. decaying states with eigenvalues corre- 
sponding to poles of the resolvent below the real energy 
axis, are thereby uncovered and can be calculated using 
standard diagonalization techniques for complex matri- 
ces. This approach is essentially exact, in the sense that 
no a priori approximations are introduced in the complex 
dilation procedure. Positions and widths of resonances 
can therefore be calculated with high precision by means 
of the complex scaling procedure [lj, Il3 • 

Quite obviously, the extension of the method to the 
Gross-Pitaevskii equation is far from straightforward 
(see, e.g., the discussion in Ref. 0]). Indeed, the pres- 
ence of the nonlinear term inhibits a direct formulation 
of complex scaling in terms of Green's functions or resol- 
vents, which, for linear systems, provide the convenient 
link to the time-dependent decay problem. Furthermore, 
self-consistent solutions of the Gross-Pitaevskii equation 
are usually defined with respect to a given normalization 
of the condensate wavefunction ip — which poses a con- 
ceptual problem for decaying states that are generally 
non-normalizable. An additional complication (which is 
quite severe from the numerical point of view, as we shall 
see later on) arises from the fact that the nonlinear term 
in the Gross-Pitaevskii contains the square modulus of 
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and is therefore a nonanalytic function of the condensate 
wavefunction. 

This last difficulty was circumvented in a recently pub- 
lished approach by Moiseyev and Cederbaum [13. In 
this work, the complex scaling transformation was ap- 
plied to the exact microscopic many-particle dynamics of 
the bosonic system, and a variational ansatz was used to 
derive from there a complex-scaled version of the Gross- 
Pitaevskii equation. This variational ansatz, however, 
was implicitly based on the assumption that the conden- 
sate wavefunction of the decaying state is entirely real, 
which effectively means that | V 7 C r ) 1 2 can be replaced by 
the analytic term tp 2 (r) in the Gross-Pitacvskii equation. 
Moiseyev and Cederbaum showed that chemical poten- 
tials and decay rates can be calculated from this com- 
plex nonlinear Schrddingcr equation, whose behaviour as 
a function of the interaction strength seems to be in qual- 
itative agreement with the intuitive expectation. 

Our ansatz in this paper is substantially different in 
the sense that we explicitly take into account the pos- 
sibility of complex-valued resonance wavef unctions. As 
pointed out above, this introduces a major complication 
of the problem, due to the presence of the |?/H r )| 2 term 
in the Gross-Pitaevskii equation. We solve this prob- 
lem by formally defining separate dilations of ip and its 
complex conjugate, and by performing explicit complex 
rotations of the wavefunction in the numerical implemen- 
tation. Our approach is supported by calculations based 
on the "real" , i.e. unsealed, Gross-Pitaevskii equation in 
presence of absorbing boundaries, with which we obtain 
good agreement. 

The paper is organized as follows: In Section III Al 
we derive the relation between complex "resonance" so- 
lutions of the stationary Gross-Pitaevskii equation and 
the actual time-dependent decay process of the conden- 
sate. For the sake of simplicity, we restrict ourselves to 
the case of the one-dimensional mean-field dynamics of a 
Bose-Einstein condensate in a magnetic or optical waveg- 
uide. The complex scaling method and its application 
to the Gross-Pitaevskii equation is explained in Section 
I11BI and details on the numerical implementation of the 
method are presented in Section Ul CI Section lTTTl contains 
the discussion of the numerical results that we obtain for 
a harmonic trapping potential with Gaussian envelopes. 



II. COMPLEX SCALING OF THE 
GROSS-PITAEVSKII EQUATION 

A. The time-dependent decay problem 

We start from the one-dimensional time-dependent 
Gross-Pitaevskii equation 

ift^U(s,t) - -^-j^*(x,t) + V(x)*(x,t) 
at 2m ox 

+2a s huj ± \^(£,t)\ 2 ^(S:,t) (1) 



which describes the mean-field dynamics of the conden- 
sate in presence of a tight cylindrical confinement with 
transverse frequency u>j_. Here, m is the mass and a s the 
s-wave scattering length of the atoms. For the sake of 
defimteness, we consider, as in Ref. |9|, the open longi- 
tudinal potential 

V(x) = -ttilu 2 x 2 exp(— ax 2 ) , (2) 

which is shown in Fig. ^ Variants of J5J were extensively 
studied in the literature on resonances (e.g. 0,0, El)- 
Since V(x) — > for x — ► ±00, the system does not 
exhibit any bound state. For not too large values of 
a, however, the condensate can be temporarily stored 
within the potential well, from where it decays via tun- 
neling through the barriers. In the linear case of non- 
interacting atoms (a s = 0), this decay process is de- 
scribed by resonance states $f(x, t) = ^(x) exp(— iEt/fi) 
where ^(x) satisifies the stationary Schrodinger equa- 
tion for the complex energy E = fi — iT/2 and exhibits 
outgoing (Siegert) boundary conditions [l!j of the form 
fy(x) — > * exp(ifc|x|) with Re(fc) > for x — > ±00 . An 
exponential decay oc exp(— Yt/fi) is therefore obtained for 
the atomic density inside the well if the system is initially 
prepared in such a resonance state. This is fundamen- 
tally different in presence of finite interaction (a s =/= 0) 
where the tunnel coupling through the barriers explicitly 
depends, via the nonlinear term in the Gross-Pitaevskii 
equation, on the density |^(ai,t)| 2 . As a consequence, 
the decay rate T varies during the time evolution of the 
condensate, which results in a nonexponential decay (see, 
e.g., 

As long as the rate T characterizing the temporal varia- 
tion of the density inside the well is comparatively small, 
the decay process can be approximately described by 
an adiabatic ansatz where the condensate is assumed 
to remain always in the energetically lowest (and most 
stable) resonance state associated with a given instan- 
taneous density |^(i,i)| 2 . To this end, we introduce 
dimensionless variables which are formally obtained by 
setting fi = m = ojq = 1. Specifically, we define the 
dimensionless potential 

v(x) — V(aox)/fiu!o — —x 2 exp(— ax 2 ) (3) 

with the dimensionless position x = x/ao and a = a^a, 
where a = y/h/ (mu ) is the oscillator length associated 
with the well (which would be of the order of several mi- 
crons for typical experimental setups). We furthermore 
define, as a function of the effective interaction strength 
g, the dimensionless resonance state ip g (x), together with 
its associated eigenvalue E g = fi g — iT g /2, as the solution 
of the stationary equation 

H(^ g (x) = E g ^ g (x) (4) 

with the dimensionless nonlinear Hamiltonian 

HW = ~-^+v(x)+g\^(x)\ 3 . (5) 
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FIG. 1: Effective longitudinal potential J3J for a = 0.1. The 
vertical dashed lines define the spatial region within which 
the wavefunction is normalized according to Eq. J7J. 



In addition to Eq. Q , we demand that tp g be normalized 
according to the condition 



m 9 



\tp g (x)\ 2 w(x)dx = 1 



(6) 



where the weight function w(x) measures the population 
inside the well. We shall use w(x) — 8(x a — x)9{x + x a ) 
in the following, i.e., 



\4>g{x)\ 2 dx. 



(J) 



where x a — 1/ \fa corresponds to the maximum of the 
barrier. We remark that the numerical results presented 
in Section ITTT1 do not sensitively depend on the particular 
choice of x a . 

A further requirement that one should impose are out- 
going boundary conditions for the wavefunction. Due to 
the selfconsistent nonlinear term in the Gross-Pitaevskii 
equation, however, this requirement cannot be formu- 
lated in the same explicit way as for the linear case 
[T^ |. Qualitatively, such outgoing boundary conditions 
imply that the local current density of ip g (x) is di- 
rected away from the trapping potential, and that no 
additional back-reflections arise outside the support of 
v(x). This means that the wavefunction evolves accord- 
ing to tpg( x ) °^ exp[i J x k(x')dx'] where the real part of 
the effective wave number k(x) is positive (negative) for 
i > i Q (x -C — x a ) and varies smoothly with x to ac- 
count for the variation of the selfconsistent potential in 
Eq. (JjJJ. We shall assume that such a condition is met 
for the wavefunction ip g . 

The adiabatic ansatz for ^(x 7 t) is now formulated as 



= y^-^g(i) (^) exp (-i J E g{tl) dt' 



(8) 



with x = x/ao and t = toot, where the effective time- 
dependent interaction strength is given by 



g(t) = 2^N(t) 
a Q luq 



(9) 



N(t) is the time-dependent population inside the well 
and decays according to the equation 

dN 



dt 



(10) 



with the initial condition iV(0) 
solved yielding 



N(t) = N exp 



No, which is formally 



T g(t') dt ' 



(11) 



As can be straightforwardly verified, ty(x,t) satisfies the 
time-dependent Gross-Pitaevskii equation as long as 
the temporal variation of ipg(uj t)( x / a o) (which is propor- 
tional to r g ) can be neglected in Eq. QJ. We thereby 
obtain 



\^(x,t)\ 2 w{x/a )dx = N(w i) . 



(12) 



The full characterization of the time-dependent decay 
process within this adiabatic picture therefore requires 
to calculate the resonance states ip g and their complex 
energies E g within the range < g < g(0) for a s > 0, and 
within g(0) < g < for a s < 0. We shall show now that 
this aim can be achieved with the method of complex 
scaling. To simplify the discussion, we shall drop the 
index g (i.e., ip = ipg and E = E g ) and consider a = 0.1 
in the following. 



B. The complex scaling transformation 

1. Linear case 

Our aim is now to calculate the resonance states for the 
nonlinear Gross-Pitaevskii equation — i.e, the solutions 
of 



with 



H(i/j)i/>(x) = Eip(x) 



H(il>) = H + g\rP(x)\ 2 , 
1 d 2 

Ho = -2dx^ +v{x) > 



(13) 

(14) 
(15) 



that exhibit outgoing boundary conditions and are nor- 
malized according to 



\tp{x)\ 2 dx = 1 



(16) 



For the linear case, i.e. in the absence of interaction (g = 
0), this aim can be straightforwardly achieved by the 
method of complex dilation. This operation is formally 
defined by the linear mapping 



4>{x) ^ i>^>(x) = R ^{x) = e w/z ^{xe w ) (17) 

which effectively amounts to a nontrivial scaling of the 
configuration space by the complex factor exp(i0). Ap- 
plying this transformation to the stationary Schrodinger 
equation fCH) (for g = 0) yields 



H^ ) ^ e \x)=E^ e \x) 



(18) 
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FIG. 2: 

for q = 
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Spectrum of the complex scaled Hamiltonian (flfjfl 
0.1. The stars represent the complex eigenvalues 



of H y "> for 9 = 0.1, which were calculated with a grid basis 
covering the range —100 < x < 100. Four resonances are fully 
uncovered at this value of the rotation angle. 



where 

Tjl@) — D TJ 7} — 1 

H — Keri<)K e — 



-3--~£ + «<-"> 0*> 



represents the (nonhermitean) complex scaled Hamilto- 
nian. 

The spectral properties of the complex Schrodinger 
equation (jT%|) are widely discussed in the literature on 
complex scaling [l(J, HH [H> 03 ■ While all eigenstates of 
the real equation 1)1 3|) will formally solve also Eq. I|18(l 
after the transformation, their normalizability properties 
might change under the operation <|17[): Bound states 
of v(x) (which are absent in our particular case) fall off 
sufficiently rapidly with increasing |x|, such that they 
remain normalizable under complex dilation (as long as 
\9\ < tt/4). This is, however, not true for continuum 
states, which turn into wavefunctions that exponentially 
diverge for |x| — * oo. For compensation, a "new" contin- 
uum of asymptotically oscillatory states arises in the ro- 
tated system, with eigenvalues along the axis E — ee~ 2l ° 
with real e. 

The rotation of the continuum axis uncovers the spec- 
tral resonances of the system, which correspond to the 
poles of the analytical continuation of the Green function 
G = (E — Ho + i<5 ) 1 to the lower half part of the com- 
plex energy plane. Those resonances turn into discrete 
complex eigenvalues E n = fx n — iT n /2 under complex 
dilation, and are represented by normalizable eigenfunc- 
tions ipffl (x) that can be straightforwardly calculated by 
diagonalizing Hq in any numerical basis. This is illus- 
trated in Fig. |21 which shows the spectrum of the rotated 
Hamiltonian i|19|) for 9 = 0.1. Besides a "continuum" of 
levels along E = ee~ 2lB (which appears as discretized due 
to the finite basis), four major resonances can be iden- 
tified below the real axis. Starting with an initial state 
that is very close to one of the resonances (i.e., which 
has, in the complex rotated system, a macroscopic over- 
lap with the associated wavefunction ?/4 ) will therefore 



lead to a decay from the well with a rate that is given 
by the negative imaginary part r„ of the corresponding 
eigenvalue. 



2. Nonlinear case 

In the nonlinear case, a major complication arises 
from the fact that the interaction-induced contribution 
g\ip(x)\ 2 to the Hamiltonian l|14|) is nonanalytic in rp. To 
avoid this complication, it is tempting to make, in a first 
approach, the replacement 



(20) 



in Eq. (|14|) . i.e. to consider the analytic nonlinear 
Schrodinger equation 



with 



H(ip)ip(x) = Etjj{x) 



H = H + g &(x)} 



(21) 



(22) 



Applying the complex scaling transformation (|17f) to 
Eq. (23} yields 



with 



(0) 



(23) 



(24) 



and H^> being defined as in Eq. I|19fl . Here the effective 
interaction strength is transformed according to 



9^ 9e = ge 



(25) 



in order to compensate the scaling prefactor e 10 / 2 in 
Eq. (|17fl . This additional scaling reflects the fact that 
g implicitly contains information about the norm of the 
wavefunction [see Eq. ©]. 

The above approach yields meaningful results (i.e. re- 
sults that are related to the actual time-dependent decay 
process of the condensate) only if the wavefunction ip(x) 
to be calculated in Eq. (0} is known to be entirely real. 
In general, this is not the case for unbound resonance 
states, which typically exhibit outgoing boundary condi- 
tions (i.e, ip( x ) exp(ifc|x|) for |x| — * oo) and are there- 
fore intrinsically complex. For such states, the replace- 
ment (|20|l represents a major modification of the prob- 
lem, and the attempt to compute the decay behaviour 
of the condensate by means of Eqs. I|23H25|I (which, in 
effect, were also used in Ref. ^5|) would lead to a predic- 
tion that is not expected to be in agreement with alter- 
native approaches based, e.g., on real-time propagation 
of the Gross-Pitaevskii equation in presence of complex 
absorbing potentials. 

To account for the fact that ip is complex-valued, 
we formally introduce a second analytic wavefunction ip 
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which coincides with the complex conjugate of ip on the 
real axis, i.e. 

ip(x) = ip*(x) for real x. (26) 

The complex scaling transformation is now defined for tp 
in the same way as for tp, 

ij>(x) i — * j} 9) (a;) = = e i6 ^(xe ie ) (27) 

with i?§ the dilation operator. The analytic continua- 
tion of the stationary Gross-Pitaevskii equation to the 
complex domain yields then 

H (e) (VO (a;) = (a) (28) 
where the complex scaled Hamiltonian is given by 

H {e) (V) = ff W + 98^ {6) (x)^ {9) (x) (29) 

with = ge~ l8 . Note that -0 is, in general, not 

identical to [tp^]*(x), the complex conjugate of ijj^ e \ It 
can, however, be obtained from ip^ via the relation 

^ {e) (x)=Re (R-e^) (x). (30) 

For weak or moderate nonlinearities (i.e., with g be- 
ing of the order of unity or smaller), the spectrum of 
self-consistent solutions of the complex scaled Gross- 
Pitaevskii Hamiltonian (|29|l can be assumend to be not 
much different from the linear case shown in Fig. [5] be- 
sides the continuum, a few quasibound resonance levels 
are expected to arise below the real energy axis. In the 
context of the actual time-dependent decay process of the 
condensate, we are mainly interested in the energetically 
lowest resonance, i.e. the one with the lowest real part /x 
of the complex eigenenergy E — /x — iT/2 (which would 
adiabatically evolve into the self-consistent ground state 
of the potential well if the barrier height was raised to 
infinity). We shall explain now how the cigcnfunction 
associated with this lowest resonance can be numerically 
calculated. 



C. Numerical calculation of the resonance state 

Our approach is based on the assumption that the low- 
est resonance state exhibits a decay rate that is not much 
larger than the decay rate of the lowest numerical con- 
tinuum states. This is indeed true in the linear case, for 
the spectrum that is shown in Fig. [21 the width of the 
lowest resonance is found to be Y/2 ~ 10~ 6 , while the en- 
ergetically lowest "continuum" state obtained from the 
numerical diagonalization (using a grid that covers the 
spatial range —100 < x < 100) decays with the rate 
r/2 ~ 10~ 4 . This observation suggests a real-time prop- 
agation approach to calculate the lowest resonance state. 
We start from a good initial approximation -0q (a;) for 



the resonant state (e.g. the harmonic eigenstate in the 
well), and numerically propagate tp^ under the time- 
dependent Gross-Pitaevskii equation 

i^\x)=H^(i, T )^\x) (31) 

in the rotated system (r is here a fictitious numerical 
"time" parameter, which is unrelated to the physical 
propagation time t in the actual decay process). In the 
above linear case, this propagation is guaranteed to con- 
verge towards the eigenstate with the lowest decay rate. 

This approach works also for the nonlinear case, but 
requires some nontrivial modifications there. On the one 
hand, the self-consistent eigenstates and eigenenergies of 
the Gross-Pitaevskii equation are defined with respect to 
a given normalization of the wavefunction. A rescaling of 
ip( 6 ^ according to the condition (|16fl is therefore required 
after each propagation step under Eq. (|31|l . On the other 

hand, ip^\x) needs to be calculated from ipf\x) in order 
to compute the nonlinear term in H^ s \ip T ). In order to 
implement the prescription (|30|) which achieves this goal, 
the complex rotation H17|) of the wavefunction needs to be 
explicitly performed. This operation, however, is known 
to be highly unstable [2lJ and requires great care in the 
numerical implementation. 

In practice, two different basis sets are simultaneously 
used in order to numerically integrate Eq. I|31|l . To per- 
form the propagation, the rotated wavefunction is ex- 
panded in a grid basis — i.e., 

™max 

4 6 \x)= £ C n (t) Xn (x) (32) 

n=— n ma x 

with 

/ n _ J 1/A* : (n - l/2)A x < x < (n + 1/2)A X 
Xn[X) \ : otherwise 

(33) 

where is a suitable grid spacing. With the standard 
finite-difference approximation for the operator of the ki- 
netic energy, 

1 9 2 Xn+l{x) +Xn-l(x) - 2Xn(x) 

"2dx^ Xn{x) - 2Al ' (34) 

we obtain a tridiagonal Hamiltonian matrix. This matrix 
is used to propagate the wavefunction according to the 
implicit scheme 

4% = (1 + if H(6) Mr)) ^ (1 - if H W (0 r ))4 e) . 

(35) 

The latter requires direct matrix-vector multiplications 
as well as the solution of linear systems of equations with 
tridiagonal matrices, which is efficiently performed by LR 
decompositions. 

For the calculation of ip , a second, nonorthogonal 
basis set is introduced in terms of the analytic Gaussian 
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functions 

^{x) = ~^=^ f-^#) (36) 

for -i/ max < v < v max . The centers x v and widths <j v 
are chosen according to 

smh(7j 

for suitable parameters 7, A x > 0, and 

0V = a!|i/+i| - x\ v \ ■ (38) 

The exponential increase of the widths a v with \v\, which 
results from Eqs. (|37[) and (|38[) . is required to ensure sta- 
ble evaluation of linear combinations of § v in the complex 
domain. ^ max is chosen such that the whole spatial range 
spanned by the grid basis is covered by the Gaussians. 

From the overlap integrals V„.„ = J Xn{x)4> v {x)dx and 
Ivy = J <f>v [x)<t>v' (x)dx, we can determine the expansion 
coefficients D— UmjBx . . . D Vmax of ip t (x) with respect to 
the Gaussian basis: 

4°\x)= J2 D u (t)M*)- (39) 

This involves again the solution of a linear system of 
equations, namely J2 v '^v,u'Du' = J2 n ^n,uC n , which is 
accomplished by the LR decomposition of the matrix 
{I v y) (the latter is effectively banded, due to the fact 
that the Gaussians fall off rapdily with increasing dis- 
tance from their center). In a very similar way, we de- 
termine the rotation of the wavefunction — i.e. the cx- 
pansion of (R-gipj. ){x) in the (unrotated) Gaussian ba- 
sis {(j)v) — by means of the overlap integrals !„.„< and 
Jv.v 1 — f >i> v {xe %e ^ v i (x)dx. Complex conjugation of the 
coefficients, the rotation Rg, and the transformation back 
to the grid basis (xn) finally lead to the state vector (C n ) 

that is associated with tp S (x) [22| . 

We should note that the above method cannot prop- 

erly reproduce ip (x) far away from the center of the 
potential. This is due to the exponential increase of 
the widths and spacings of the Gaussians, which pre- 
vents a perfect representation of rapidly oscillating wave- 
functions far away from the origin. We remark, how- 
ever, that in the time-dependent Gross-Pitaevskii equa- 

tion 131|) tp (x) is multiplied by the square of ip (e) (x), 

and the term {tp if>^ip^)(x) decreases exponentially 
with \x\. Hence, apart from the special case of very large 
g and rapidly decaying resonance states, we obtain an 

approximation for tp which is reasonably good in the 
relevant spatial regime to ensure stable convergence of 
our algorithm. 
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FIG. 3: Chemical potentials (upper panels) and decay rates 
(lower panels) of the lowest resonance state with even and 
odd parity (left and right column, respectively), as a func- 
tion of the effective interaction strength g. The inserts show 
the spatial density of the corresponding resonance states at 
g = 1 (with the trapping potential indicated by dashed lines). 
Converged calculations of /j, and F can be performed up to 
interaction strengths g where the chemical potential reaches 
the maximum barrier height of the potential, marked by the 
horizontal dashed lines in the upper panels. For g < —1.1, 
the lowest even state becomes stable, which manifests itself 
in a negative chemical potential and a vanishing decay rate. 



III. NUMERICAL RESULTS 
A. Chemical potentials and decay rates 

Fig-Elshows the chemical potentials and decay rates of 
the energetically lowest and of the first excited (soliton- 
like) quasibound state as a function of g (lower and up- 
per solid curves, respectively) for a = 0.1. Due to the 
symmetry of the potential, both states can be calculated 
with the implicit propagation scheme described above, 
where the parity of the state is selected by the choice of 
the initial wavefunction. Typical parameters employed 
in the calculations are A x ~ 0.01 for the grid spacing, 
— 100 < x < 100 for the spatial range covered by the grid, 
9 = 0.05 for the complex scaling angle, and 7 ~ 0.05 for 
the parameter characterizing the exponential increase of 
the width of the Gaussians. 

In the regime of small and moderate nonlincarities 
(Is I S 1), the numerical procedure described in the pre- 
vious subsection is further optimized: To avoid time- 
consuming rotations of wavefunctions as best as possible, 
the renormalization of if)™' as well as the calculation of 
the nonlinear term in the Gross-Pitaevskii Hamiltonian 
are performed only every 100th propagation step accord- 
ing to Eq. |J5SJ. We verified, in any case, the convergence 
of the calculation by a criterion that is independent of the 
real-time propagation approach, namely by the require- 
ment that the norm of 8^ = - fall 
below a given precision limit. To obtain highly accurate 
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MCS /^ABC MCAP 



0.4601 
0.7954 
1.0765 
1.3190 
1.5315 
1.7236 
1.9043 



0.4602 
0.7955 
1.0773 
1.3193 
1.5313 
1.7230 
1.9034 



0.4602 
0.7954 
1.0772 
1.3192 
1.5312 
1.7231 
1.9035 



Tcs/2 Tabc/2 Tcap/2 



9.35e-7 9.55e-7 9.62e-7 

1.82e-5 1.81e-5 1.80e-5 

1.55e-4 1.56e-4 1.56e-4 

8.05e-4 8.04e-4 8.05e-4 

2.76e-3 2.75e-3 2.75e-3 

6.65e-3 6.66e-3 6.63e-3 

1.24e-2 1.23e-2 1.23e-2 



TABLE I: Chemical potentials and decay rates of the low- 
est quasibound state, calculated with the complex scaling ap- 
proach (CS) and with real-time propagation methods using 
absorbing boundary conditions (ABC) and complex absorb- 
ing potentials (CAP). 



results for the chemical potentials and the decay rates T, 
it is convenient to perform, for a given value of g, sev- 
eral calculations with different numerical grid spacings 
A x (where the total extent of the grid is kept constant), 
and to extrapolate /i and T in the limit — > 0. In this 
way, rather low decay rates down to T ~ 10~ 7 can be 
reproduced. 

Since no assumption about the sign of g was made 
in our approach, we can also calculate resonance (and 
bound) states for the case of attractive interaction. As 
we see in Fig. |3 the lowest even state of the trapping po- 
tential undergoes, for increasing attraction between the 
atoms, a transition from a resonance to a bound state, 
which is manifested by a negative chemical potential and 
a vanishing decay rate. In agreement with the calculation 
by Moiseyev et al. , which were performed for the same 
potential © with a = 0.2, we find that this transition 
occurs at g ~ —1.1. 

For strong nonlinearities (\g\ > 1), the renormaliza- 
tion of tfi^ and the update of the nonlinear term need 
to be performed after each propagation step in order 
to ensure stable convergence of our approach. A fur- 
ther modification is introduced in the regime of large de- 
cay rates V ~ 0.01, i.e. where fi approaches the height 
of the potential barriers: To avoid unwanted conver- 
gence to energetically low continuum states (with de- 
cay rates T ~ 10~ 4 , see Sec. Ill Bp . we replace, in this 
regime, the real-time propagation approach by a prop- 
agation along a complex time path, given by r = fe lip 
with real f and < <p < 26. A stationary solution 
tp of the Gross-Pitaevskii equation will then evolve ac- 
cording to \ipr\x)\ = exp[(iJ,tp — r/2)f] \ipQ e \x)\ for 
small tp, and quasibound states with finite /i will thereby 
be more strongly enhanced than low-energy continuum 
states close to the threshold. Indeed, we find that this 
modification allows converged calculations of resonances 
with chemical potentials close to the barrier height. 

Table IJ shows a comparison of the chemical potentials 
and decay rates with the ones that are obtained from real- 
time propagation approaches based on the original (i.e., 
unrotated) Gross-Pitaevskii equation In those calcu- 
lations, the decay through the boundaries of the numer- 




FIG. 4: (Color online) Wavefunctions of the even and odd res- 
onance states at g = 1 (upper and lower panels, respectively; 
solid lines: real part; dashed lines: imaginary part). The 
wavefunctions were calculated for the parameters 7 = 0.05 
and A x = 0.1 of the Gaussian basis set (see Eq. 1371 '). The 
magnifications in the right panels show that the outgoing tails 
of the wavefunctions are faithfully reproduced till x ~ 30. 



ical grid is accounted for by absorbing boundary condi- 
tions [23L |24| as well as by a complex absorbing potential 
of the form T4bs(^) = — i(x/100) 10 (for the spatial range 
— 100 < x < 100). As in the complex scaling approach, 
the wavefunction was propagated by an implicit scheme 
like (|35f) (now under the real, unsealed Gross-Pitaevskii 
Hamiltonian) and renormalized after each step to sat- 
isfy the condition J7J. The agreement between the three 
methods is reasonably good. 

We should note, however, that strongly deviating val- 
ues for r would generally be obtained if the resonance 
states were calculated according Eqs. (|23I24[) under the 
assumption that the wavefunction -0 is entirely real. For 
g = 1, e.g., the real-time propagation approach (|35|l us- 
ing [from Eq. $ZQ ] instead of yields essentially 
the same chemical potential fi ~ 0.795, but a different 
decay rate T/2 ~ 5 x 10~ 5 . Also for g = 5 we obtain 
r/2 ~ 1.6 x 10~ 2 which considerably overestimates the 
actual rate r/2 ~ 6.7 x 10~ 3 . The above calculations 
based on were performed under the slightly modi- 
fied side condition J"^™ [tp(x)] 2 dx — 1, which would be 
consistent with Eq. I|16|) for purely real ip [the propaga- 
tion of tp^ under does not seem to converge with 
the original side condition (|16l) ]. 



B. Wavefunctions 

While the chemical potentials fj, and decay rates T are 
apparently well reproduced by the method of complex 
scaling, the wavefunctions ip(x) of the resonance states 
(defined according to Eq. (|13H . i.e., in the unrotated 
framework) can be calculated only in a limited spatial 
regime in the vicinity of the potential well. We attribute 
this to the exponential increase of the widths u v and cen- 
ters x v with \v\, which was introduced in order to ensure 
stable back-rotation to the real domain. As already men- 
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FIG. 6: Time-dependent decay of the condensate. Plotted is 
the number N of atoms inside the potential well as a func- 
tion of time, with No = N(t — 0) such that g — 1 initially. 
The dashed line shows the corresponding decay behaviour of 
noninteracting atoms. 



FIG. 5: (Color online) Wavefunction of the even resonance 
state at g = 5 (solid line: real part; dashed line: imagi- 
nary part), calculated with the complex scaling method (up- 
per panel) and with a real-time propagation approach using 
complex absorbing potentials (lower panel). The comparison 
shows that the complex scaling approach fails to correctly re- 
produce the resonance wavefunction for \x\ > 15, even though 
the chemical potential and the decay rate are in agreement for 
both methods (see TableHJ. The insert in the lower panel dis- 
plays the density of the resonance state continued till x = 500, 
the tail of which was computed by integrating the free one- 
dimensional Gross-Pitaevskii equation (with v(x) = 0) with 
the complex effective chemical potential [i g=5 — iT g= $/2, start- 
ing from i/'CAp(a; = 10). A singularity of \ip\ 2 is encountered 
at x ~ 485 where the density exceeds the critical value fi/g. 



tioned, this exponential scaling prevents a perfect repre- 
sentation of rapidly oscillating wavefunctions far away 
from the origin. For low and moderate nonlinearities 
(.9 S 1), we typically reproduce ip(x) up to distances of 
the order of |x| ~ 50 from the center of the potential. 
This is illustrated in Fig. 0] where the wavefunctions of 
the even and odd resonance states at g = 1 are plotted. 

For large g, the spatial range in which ip can be cor- 
rectly reproduced becomes shorter. This is demonstrated 
in the upper panel of Fig. |3]which displays the wavefunc- 
tion tpcs of the even resonance state at g = 5 calculated 
with the complex scaling approach. In contrast to the ac- 
tual resonance wavefunction (shown in the lower panel of 
Fig. [5J, the density \4>cs(x)\ 2 exhibits unphysical max- 
ima near |x| ~ 25 and decreases rapidly for larger \x\. 
We believe that those features arise due to the imperfect 

computation of the complex conjugate function ip (x) 
by means of explicit rotations within the Gaussian basis 
set. It should be noted, however, that the values for \x 
and r/2 that are extracted from this wavefunction are 
nevertheless in good agreement with the chemical poten- 
tials and decay rates calculated with alternative methods 
at g = 5 (see Table HJ. 

Let us finally remark that from a formal point of 
view, even the "true" resonance wavefunction cannot 
be defined within an infinitely large spatial regime for 



g > 0: Due to the outgoing boundary conditions, the 
density ^(a;)! 2 of the resonance wavefunction constantly 
increases with increasing |x| outside the potential well, 
until it exceeds the critical value fj,/g where the kinteic 
energy formally vanishes. Beyond that point, a rather 
sudden increase of the density towards infinity is en- 
countered [25|, which indicates that permanently time- 
dependent behaviour would be expected for the actual 
decay process in that spatial regime. We should note, 
however, that the associated critical distance |x| = x c at 
which |i/)(a;)| 2 equals [if g is much larger than the spatial 
regime taken into account in our numerical calculations; 
for the even resonance state at g = 5, e.g., we obtain 
x c ~ 500 from fj, ~ 1.7 and r/2 ~ 0.007, while even 
larger critical distances (x c ~ 10 6 ) would be expected for 
the more long-lived resonance states around g = 1. 



C. The time-dependent decay process 

With the above method to calculate the decay rates 
of the resonance states for given nonlinearity strengths 
g, we are now in a position to predict the actual, real- 
time decay process of the condensate. As was derived 
in Section [II Al from the adiabatic approximation (jSJl, the 
number of atoms in the potential well decreases with time 
according to 

j t N(t) = -T g(t) N(t) , (40) 

where g(t) is proportional to N(t) via Eq. |J3J|: 

g(t)=2^^N(t). (41) 

Multiplying Eq. I|40|) with the constant prefactors appear- 
ing in Eq. (|41|l yields the ordinary differential equation 



dt 



g(t) 



g(t) 



git). 



(42) 



This equation can be efficiently integrated e.g. with a 
Runge-Kutta solver, provided the rates T g are known in 
the range < g < g(0). 
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Fig. [U shows the result of this integration for the case of 
repulsive interactions (a s > 0). Here, the initial number 
of atoms N(Q) was assumed such that g(0) = 1. The 
calculation was based on the decay rates To, To.i . . .T\ 
that were computed with the complex scaling method, 
and employed cubic interpolation to obtain intermediate 
values for T g . We clearly see that N(i) decays initially 
according to a subexponential law, due to the fact that 
r decreases with decreasing g (a superexponcntial law 
would be encountered for attractive interaction). After 
about 10 5 time units, the exponential decay behaviour 
of a noninteracting condensate (with decay rate Tq ~ 
2 x 10 ~ 6 ) is recovered. 

For comparison, we also performed a full time- 
dependent calculation of the decay process of the conden- 
sate. To this end, we first computed, using again Eq. 
as side condition, the self-consistent ground state of the 
condensate at g = 1 in the "closed" trapping potential 



[ v(x a ) : \x\ > x a 

with x a = 1/y/a, where the exterior part of v{x) is re- 
placed by a constant level corresponding to the height 
of the barriers. At t = 0, v(x) was suddenly "opened", 
i.e. replaced by v(x), and the wavefunction of the con- 
densate's ground state was propagated under the time- 
dependent Gross-Pitaevskii equation at g = 1, where 
absorbing boundary conditions were used to account 
for the decay through the boundaries of the numerical 
grid. The agreement between the two approaches is ex- 
tremely good. At t — 4 x 10 5 for instance, when the 
surviving population of the condensate has decayed to 
N(t) ~ O.lAo, we find that the relative difference be- 
tween the above propagation method and the integration 
of Eq. I|42|l (which takes much less CPU time, including 
the calculation of the decay rates T g ) is of the order of 
AN/N(t) ~ 0.01. This shows that the rates T g obtained 
from the complex scaling approach are sufficiently pre- 
cise to predict the time-dependent decay dynamics of the 
condensate. 



IV. CONCLUSION 

In summary, we have shown that the method of com- 
plex scaling can be used to calculate metastable mean- 



field states of Bose-Einstein condensates that are con- 
fined in trapping potentials with finite tunneling bar- 
riers. Our approach is based on the complex rotation 
x xe l6 of the Gross-Pitaevskii equation and employs, 
as key ingredients, separate dilations of the condensate 
wavefunction ip and its conjugate %p, together with an ad- 
ditional scaling g i— > ge~ %e of the interaction strength. A 
real-time propagation approach is used to calculate the 
lowest symmetric and antisymmetric quasibound states 
of the trapping potential and to determine the associ- 
ated chemical potentials and decay rates. We find good 
agreement with alternative propagation methods that are 
based on the unrotated Gross-Pitaevskii equation and 
use absorbing boundary conditions as well as complex 
absorbing potentials to describe the decay. 

Though only exemplified for the harmonic trap with 
Gaussian cutoffs, our method is sufficiently general to be 
applied to other trapping geometries, such as the dou- 
ble barrier potentials that are proposed in the context 
of resonant transport [2~j. |2(| . This is indeed confirmed 
in a recent study on the nonlinear Wannier-Stark prob- 
lem, where the decay rates of the condensates were repro- 
duced with the complex scaling method |27[ • The formal- 
ism presented here can, in addition, be straightforwardly 
generalized to describe decay processes in two- and three- 
dimensional geometries (which would naturally involve 
a higher numerical effort). Moreover, the framework of 
complex scaling seems also suited to be applied beyond 
the pure mean-field description of the condensate, e.g. in 
the context of fragmentation |28| and within the micro- 
scopic quantum dynamics approach j2^. We therefore 
expect that the complex scaling technique might develop 
into a useful tool to predict the storage time of conden- 
sates in microscopic trapping potentials, and to address 
the issue of resonances of the nonlinear Gross-Pitaevskii 
equation from the conceptual point of view. 
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